\section{Estimate $S_t$ by a L\'evy process}
In earthquake modeling, one lets $S_T$ be the "stress potential" at time $t$ and $s$ be the "threshold value". At time $\tau = \inf_{t}\{S_T>s\}$ an earthquake occurs, and the potential is reset to zero.\\

\noindent
In between earthquakes the stress potential builds up according to a pure jump L\'evy process with L\'evy measure on $[0,\infty)$ given by

$$
d\nu(x)=\left(\frac{1}{1+x}\right)^3dx, \quad x \ge 1.
$$

\noindent
In this problem we assume that the threshold value is $s=20$.

\subsection{Produce a histogram of the severity of an earthquake}
We want to produce a histogram of the severity of an earthquake, which means that our task is to simulate the amount with which $S_T$ exceeds $s=20$ at the time $\tau$.\\

\noindent
In order to simulate $S_T$ we need to know something about the L\'evy process. We begin by showing that $\nu$ has finite mass.

$$
\int_{1}^{\infty} (1+x)^{-3} dx = \brackk{\tfrac{-1}{2}(1+x)^{-2}}_1^\infty = \frac{1}{8}
$$

\noindent
Then we know that $S_T$ is a compound Poisson process 

$$
S_t = \sum_{i=1}^{N\paren{t}}\xi_i,
$$

\noindent
where $N(t)$ is a Poisson process with $\lambda=\tfrac{1}{8}$, and the jump sizes, $\xi_i$, are iid. random variables with distribution $\frac{\nu(dx)}{\lambda} = 8(1+x)^{-3}dx$ on the interval $[1,\infty)$.\\

\noindent
We want to use the inverse transform method to simulate jump sizes, and therefore we need to know the distriution function
$$
F(y) = \int_1^y8(1+x)^{-3}dx = 1 - 4(1+y)^{-2}.
$$
\noindent
From this the inverse is easily found to be $F^{-1}(u) = 2 (1-u)^{-\tfrac{1}{2}}-1$.\\

\noindent
We notice that since $S_T$ is a compound Poisson process the time at which it exceeds the threshold must be a jump time, and since we are only interested in the value with which it exceeds the threshold we need not simulate the actual jump times but only the jump sizes. Then we add up the jump sizes and the result is the difference between $S_t$ and $s$ the first time $S_t$ exceeds $s$.\\

\noindent
This method is summarized in algorithm~\vref{Problem4Algo} and in appendix \ref{code4a} there is an implementation in \texttt{R} that generates 1000 samples and plots the following histogram.

\begin{center}
\includegraphics[width=10cm]{problem4_1.eps}	   
\end{center}

\begin{algorithm}
\caption{Generate values from $S(\tau)-s$}
\label{Problem4Algo}
\begin{algorithmic}[1]
  \STATE $s \leftarrow 20$
  \STATE $S \leftarrow 0$
  \STATE
  \REPEAT
  \STATE Generate $u$ $\sim$ Unif(0, 1)
  \STATE $S \leftarrow S + 2 (1-u)^{-\tfrac{1}{2}}-1$
  \UNTIL{$S > s$}
  \STATE 
  \RETURN $S - s$
\end{algorithmic}
\end{algorithm}

\subsection{European call option}
Here we consider the standard European call option, where the objective is to estimate

$$
E[e^{-rT}(S_T-K)^+],
$$

\noindent
where $S_T$ denotes the stock price at time $t$. We assume that we model 

$$
\log S_t - \log S_0 = \left(r-\frac{\sigma^2}{2} \right)t+\sigma W_t + J_t,
$$

\noindent
where $W_t$ is a Brownian motion, and $J_t$ is a pure jump process with L\'evy measure 

$$
d\nu(x) = \frac{1}{5}\frac{x^{\tfrac{1}{2}}+3}{x^{\tfrac{7}{4}}}dx, \quad x>0.
$$

\noindent
First we look at the pure jump process and see that $\nu$ has infinite mass in this case.

\begin{eqnarray*}
\int_0^\infty d\nu(x) &=& \frac{1}{5}\int_0^\infty \frac{x^{\tfrac{1}{2}}+3}{x^{\tfrac{7}{4}}}dx\\
 &=& \frac{1}{5}\int_0^\infty x^{-\tfrac{5}{4}}+3x^{-\tfrac{7}{4}}dx\\
 &=& \frac{1}{5}\left[-4x^{-\tfrac{1}{4}}-3\tfrac{4}{3}x^{-\tfrac{3}{4}}\right]_0^\infty\\
 &=& \frac{-4}{5}\left[x^{-\tfrac{1}{4}}+x^{-\tfrac{3}{4}}\right]_0^\infty\\
 &=& \infty
\end{eqnarray*}

\noindent
This means that we can't do the straight forward compound Poisson simulation. Instead we divide $J_T=J^{(1)}(t)+J^{(2)}(t)$ where

\begin{eqnarray*}
d\nu^{(1)}(x) &=& 1_{x<1}d\nu(x) \\
d\nu^{(2)}(x) &=& 1_{x\ge 1}d\nu(x)
\end{eqnarray*}

Now its's easily shown that $J^{(2)}(t)$ has finite mass

\begin{eqnarray*}
\int_0^\infty d\nu^{(2)}(x) &=& \int_1^\infty d\nu(x)\\
 &=& \frac{-4}{5}\left[x^{-\tfrac{1}{4}}+x^{-\tfrac{3}{4}}\right]_1^\infty\\
 &=& \frac{8}{5}
\end{eqnarray*}

such that it can be simulated by a compound Poisson process 

$$
J^{(2)}(t) = \sum_{i=1}^{N(t)}\xi_i,
$$

\noindent
where $N(t)$ is a Poisson process with $\lambda=\tfrac{8}{5}$, and the jump sizes, $\xi_i$, are iid. random variables with distribution $\frac{\nu^{(2)}(dx)}{\lambda}$ on the interval $[1,\infty)$.\\

\noindent
We want to use the inverse transform method to simulate jump sizes, and therefore we need to know what the distriution function

\begin{eqnarray*}
F(y) &=& \frac{5}{8}\int_1^y\frac{1}{5}\frac{x^{\tfrac{1}{2}}+3}{x^{\tfrac{7}{4}}}dx\\
 &=& \frac{5}{8}\frac{-4}{5}\left[x^{-\tfrac{1}{4}}+x^{-\tfrac{3}{4}}\right]_1^y\\
 &=& \frac{-1}{2}\left(y^{-\tfrac{1}{4}}+y^{-\tfrac{3}{4}}-2\right)\\
 &=& 1 - \frac{1}{2}\left(y^{-\tfrac{1}{4}}+y^{-\tfrac{3}{4}}\right)
\end{eqnarray*}

\noindent
We don't want to find the inverse function, but it is easily seen that $F(y)$ is strictly increasing on the interval $[1, \infty)$. This means that the inverse function exists, and for an $u$ $\sim$ Unif(0,1) we can find it numerically.\\

\noindent
In order to simulate $J^{(1)}(t)$ we calculate 

\begin{eqnarray*}
\mu_1 &=& \int_0^1 x d\nu(x)\\
 &=& \frac{1}{5}\int_0^\infty x^{-\tfrac{1}{4}}+3x^{-\tfrac{3}{4}}dx\\
 &=& \frac{1}{5}\left[\tfrac{3}{4}x^{\tfrac{3}{4}}+3\tfrac{1}{4}x^{\tfrac{1}{4}}\right]_0^1\\
 &=& \frac{3}{10}\\
\sigma_1^2 &=& \int_0^1 x^2 d\nu(x)\\
 &=& \frac{1}{5}\int_0^\infty x^{\tfrac{3}{4}}+3x^{\tfrac{1}{4}}dx\\
 &=& \frac{1}{5}\left[\tfrac{7}{4}x^{\tfrac{7}{4}}+3\tfrac{5}{4}x^{\tfrac{5}{4}}\right]_0^1\\
 &=& \frac{3}{5}
\end{eqnarray*}

\noindent
We can now simulate $J^{(1)}(t)$ by using the property that 

$$
\frac{J^{(1)}(t)-\mu_1t}{\sigma_1} \overset{d}{=} Z_t
$$

\noindent
where $Z_t\sim$ N(0,t).\\

\noindent
This can now be combined into a procedure that estimates the expectation of the European call option. This is done by first noticing that the only value of $S_t$ that we need is for $t=T$. The we see that this value can be calculated as 

$$
S_T=S_0e^{\log S_T-\log S_0}=S_0e^{\left(r-\frac{\sigma^2}{2} \right)T+\sigma W_T + J_T},
$$

\noindent
where $W_T=W_T-W_0\sim$ N(0, T), and $J_T=J_T^{(1)}+J_T^{(2)}$ can be estimated by the procedure described above. Now we combine this into algorithm \ref{Problem4bAlgo}. This algorithm is implemented in the code in appendix \ref{code4b}.\\

\begin{algorithm}
\caption{Estimate the mean of the standard European call option}
\label{Problem4bAlgo}
\begin{algorithmic}[1]
%  \STATE $r \leftarrow 0.04$
%  \STATE $\sigma \leftarrow 0.3$
%  \STATE $S_0 \leftarrow 100$
%  \STATE $T \leftarrow 1$
%  \STATE $K \leftarrow 105$
%  \STATE 
%  \STATE $N \leftarrow 10000$
  \STATE $\hat{\mu}\leftarrow 0$
  \FOR {$n \leftarrow 1$ to $N$}
  \STATE $j^{(2)} \leftarrow 0$
  \STATE $t \leftarrow 0$
  \REPEAT
  \STATE Generate $t_i$ $\sim$ Exp($\tfrac{8}{5}$)
  \STATE $t \leftarrow t + t_i$
  \IF{$t \le T$}
  \STATE Generate $u$ $\sim$ Unif(0,1)
  \STATE Find $\xi_i$ by root search for $y\ge 1$ such that $0=1 - \frac{1}{2}\left(y^{\tfrac{-1}{4}}+y^{\tfrac{-3}{4}}\right)-u$
  \STATE $j^{(2)} \leftarrow j^{(2)} + \xi_i$
  \ENDIF
  \UNTIL{$t \ge T$}
  \STATE 
  \STATE Generate $x$ $\sim$ Norm(0, T)
  \STATE $j^{(1)} \leftarrow (x + \frac{3}{10}T)\sqrt(\frac{3}{5})$
  \STATE 
  \STATE Generate $w$ $\sim$ Norm(0, T)
  \STATE 
  \STATE $S_T \leftarrow S_0e^{\left(r-\frac{\sigma^2}{2} \right)T+\sigma w + j^{(1)} + j^{(2)}}$
  \STATE $\hat{\mu}\leftarrow \hat{\mu} + \max(S_T-K,0)$
  \ENDFOR
  \STATE $\hat{\mu} \leftarrow \hat{\mu}e^{-rT}/N$
  \RETURN $\hat{\mu}$
\end{algorithmic}
\end{algorithm}

\noindent
This algorithm doesn't work because the simulated $J_T^{(2)}$'s become very large. This is because the distribution function in the simulation of the jump sizes in $J_t^{(2)}$ increases very slowly. This can be seen in the graph below where we have plotted $y$ against the distribution function $F(y)$ for values of $y$ between 1 and 1000000. We have also added lines for $F(y)=0.99$ and $F(y)=0.98$ to show that more than 1 percent of the simulated $y$'s will be larger than 1000000.
 
\begin{center}
\includegraphics[width=12cm]{problem4_2.eps}	   
\end{center}

\noindent
Jump sizes of up to and over 1000000 in $J^{(2)}(T)$ makes our algorithm numerically unstable, which means that we can't produce the estimate that we set out to produce.

\newpage
